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Abstract 


The membrane electrode assembly (MEA) pressure distribution is an important factor that affects the performance of polymer membrane 
electrolyte fuel cell (PEMFC) stack. However, the general rules for assembly parameters that affect the MEA pressure distribution are hardly 
reported. In this study, a robust design analysis based on response surface methodology (RSM) was performed on a simplified fuel cell stack in 
order to identify the effect of assembly parameters on the MEA pressure distribution. The assembly pressure and bolt position were considered 
as randomly varying parameters with given probabilistic property and acted as the design variables. The max normal stress and normal stress 
uniformity of the MEA were determined in terms of the probabilistic design variables. The reliability of the robust design has been verified by 
comparing the robust solution with the optimal solution and an arbitrary solution. 


© 2007 Elsevier B.V. All rights reserved. 
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1. Introduction 


Fuel cells are recognized as potentially environmentally 
friendly power sources for residential, portable and transporta- 
tion applications [1]. The interest in the usage of PEMFC in a 
variety of power generation applications is increasing every day, 
which makes the performance improvement, economical and 
energy efficiency of fuel cell necessary. In a single PEMFC, a 
MEA is sandwiched between two bipolar plates housing the flow 
channels. Multiple fuel cells are stacked together in most prac- 
tical applications to provide sufficiently high power and desired 
voltage. This configuration results in the amplification of the 
losses from contact resistance between contacting components 
[B]. 

According to the former researches, one of the most impor- 
tant factors that affect the contact resistance is the component 
pressure distribution of PEMFC stack. For PEMFC stacks 
using common graphite bipolar plates which are not flexible, 
increasing the pressure on MEA leads to increasing the electric 


* Corresponding author. Tel.: +86 2134206303; fax: +86 2134206340. 
E-mail address: xmlai @sjtu.edu.cn (X. Lai). 


0378-7753/$ — see front matter © 2007 Elsevier B.V. All rights reserved. 
doi: 10.1016/j.jpowsour.2007.05.066 


conductivity and reducing the permeability of the assembly [3]. 
However, the brittle gas diffusion layer can be damaged if too 
much pressure is applied. There have been some PEMFC stacks 
reported using flexible bipolar plates. For example, Yan et al. [4] 
developed a type of cheap expanded graphite plate material and 
a production process for fuel cell bipolar plates. They success- 
fully assembled 1 and 10 kW stacks using the expanded graphite 
bipolar plates. Hwang and Hwang [5] assembled a double-cell 
PEMFC using Grafoil™ flow-field plates and studied the perfor- 
mance of the stack. For these flexible bipolar plate designs, not 
only the MEA pressure distribution but also the bipolar plates 
pressure distribution are very important. To some extent, the 
latter is more important. In this study, we mainly focus on the 
PEMFC stacks using common graphite plates. As to these stacks, 
because of the relatively thin dimensions and low mechanical 
strength of the MEA versus bipolar plates and end plates, and 
the requirement of small contact resistance, the most important 
goal in the stacking design and assembly is to achieve a proper 
and uniform MEA pressure distribution [6]. 

Lee et al. [6] proposed a finite element analysis (FEA) model 
and analyzed the MEA pressure distribution under given assem- 
bly pressure. They demonstrated that the point-load stacking 
design was not a good method for obtaining a uniform pressure 
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distribution. Lee et al. [7] reported the changes in fuel cell per- 
formance as a function of the compression pressure resulting 
from torque on the bolts that clamped the fuel cell. Yoon et al. 
[8] and Mishra et al. [2] conducted some experiments on mea- 
suring the contact resistance under different assembly pressure. 
Vlahinos et al. [3] first gave a model considering the material 
and manufacturing variations. They analyzed the effects of these 
variations on MEA pressure distribution. Zhou et al. [9] stud- 
ied the influence of the clamping force on the performance of 
PEMFC and made the conclusion that there would be a maxi- 
mum power density if an optimal clamping force is found for 
a practical fuel cell system. For other fuel cell types, such as 
the solid oxide fuel cell stacks, the pressure distribution also has 
effect on the contact resistance. Koch and Hendriksen [10] inves- 
tigated the load behavior of the contact resistance and found a 
power law dependence between the load and the contact resis- 
tance on a solid oxide fuel cell stack. The conclusion was drawn 
that contact resistance over an interface is highly dependent on 
the contact load. 

Nevertheless, the general rules for assembly parameters that 
affect the MEA pressure distribution have been seen in almost 
none report. For example, the amount of assembly pressure has 
been always determined by the trial-and-error process. On the 
other hand, the end plate assembly bolts may be located at the 
four corners or at the middle of four margins of the end plate. 
Furthermore, in reality, there are always variations in the assem- 
bly parameters. It is neither physically possible nor financially 
feasible to completely eliminate the variations. Therefore, the 
need exists to identify the robust solution for assembly parame- 
ters that can produce the best MEA pressure distribution and are 
also insensitive to the variations in the assembly parameters. 

In this study, it is aimed to investigate the effects of the assem- 
bly pressure and the position of end plate bolts on the MEA 
pressure distribution. A numerical model of a PEMFC stack was 
developed and the robust solution of the assembly parameters 
was achieved through the methodology described in this paper. 


Numerical modeling procedure 


CAD model 


Material properties 


Boundary conditions 
Loading conditions 
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| 
Model assumptions | 


Establish FEM model 


2. Methodology 


The main purpose of this research is to develop a methodol- 
ogy and a FEA simulation procedure to establish the numerical 
tools for the evaluation of the stack design and cell assembly 
parameters. The schematic plot of the methodology is shown in 
Fig. 1, which is composed of a numerical modeling procedure 
and a robust design process. 

The well-established finite element method was employed 
for the numerical model. At first, the dimensions of all fuel cell 
components were collected to construct the CAD model. Then, 
the mechanical properties, the loading and boundary conditions, 
and the behavior of the contacting interface of the components 
must be consistent with the actual physical situation. Finally, the 
proper types of elements for each component and their interfaces 
must be selected to allow a realistic physical behavior. Meshing 
is also important to obtain an accurate result. The significant dif- 
ference in thickness between the components requires a special 
consideration in the meshing scheme. During the creation of the 
finite element method (FEM) model, all the dimensions, load- 
ing and boundary conditions, and mechanical properties are fully 
parameterized. Since the type and amount of assembly pressure 
depends on the type of stacking designs, the FEM model devel- 
oped above is designed to be able to simulate different stacking 
designs and compression methods. For example, by applying the 
assembly pressure on the top surface of the end plate, the model 
can be used to simulate the “surface-load” stacking design. By 
adding the number of assembly bolts, the model can also be used 
to simulate the “line-load” stacking design. 

In order to find the robust solution of the assembly parame- 
ters, a robust design process was developed. The RSM was used 
to establish the response variables in a form of second-order 
polynomial of the design variables. After establishment of the 
second-order polynomial, the transmitted variation in the design 
variables can be obtained using the error propagation equation 
(POE). At last, an overall desirability function was formed by 
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Fig. 1. Schematic plot of the simulation methodology. 
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combining the single desirability function of each response. By 
maximizing the overall desirability function, the robust solution 
was obtained. After verification of robustness, the established 
numerical simulation procedures can be used to evaluate a new 
stacking design and/or optimize cell assembly parameters. 


3. Numerical model for a PEMFC stack 


A typical PEMFC stack consists of a pair of end plates, sev- 
eral bipolar plates and several MEAs. Each MEA is sandwiched 
between two bipolar plates. On the end plate, four bolts hold the 
stack together by imposing certain pressure. 

In this study, a two-cell stack model was developed and the 
FEM model is shown in Fig. 2. The commercial code of ANSYS 
was used to build the FEM model. The 3D elements solid 45 and 
solid 95 were used to represent the bipolar plate, the MEA and 
the end plate, where solid 45 and solid 95 are popular element 
types that can be used for 3D modeling of solid structures. Solid 
45 is defined by eight nodes having three degrees of freedom 
at each node and solid 95 defined by 20 nodes having three 
degrees of freedom per node is a higher-order version of solid 
45. A combination of mapped meshing and automatic meshing 
was adopted in order to ensure proper element connectivity and 
a correct aspect radio. The interfacial nodes between the bipolar 
plate and the MEA, and between the end plate and the bipolar 
plate were coupled in the normal direction to model the contact 
behavior as shown in Fig. 2(A). 

The assembly pressure was applied through four bolt and 
nut assemblies. To simplify the numerical model, the bolts and 
nuts were ignored in the finite element model. The assembly 
pressure was applied directly at the contacting areas of the end 
plate as shown in Fig. 2(B). The FEM model has to be properly 
constrained in order to prevent the free movement. Due to the 
symmetry of the stack, half of the stack was employed in the 
analysis and the symmetrical boundary conditions were applied 
at the bottom side of the model. 

The bolt position and the assembly pressure are both param- 
eterized. The elements are color-coded based on their material 
properties. The dimension and mechanical property of each com- 
ponent is listed in Table 1. The mechanical properties of the 
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Fig. 2. (A and B) Half of the finite element model of a two-cell PEMFC stack. 


Table 1 
Dimension and mechanical property of each component [3] 


Component Modules of Poisson’s Size (mm) 
elasticity (MPa) ratio 

End plate 70,000 0.3 84 x 84 x 8 

MEA 21 0.001 50 x 50 x 0.457 

Bipolar plate 5,100 0.3 50 x 50 x 1.27 


components are cited from Ref. [3]. In practice, bipolar plates 
and MEA can be anisotropic due to either their structures or man- 
ufacturing processes. However, our simulation is a 3D structure 
analysis and the anisotropic material behavior is complicated 
and time consuming for the simulation. Thus, to simplify the 
model, the material behavior was assumed isotropic. The model 
involved the following assumptions: 


(1) The small fillets were ignored. 

(2) The effect of gravity was neglected. 

(3) The material behavior was assumed to be linear elastic and 
isotropic. 


4. Robust design process for assembly parameters 
4.1. Design and response variables definition 


As shown in Fig. 3, the bolt position can be expressed as 
the vertical position dimension Pos, and the horizontal posi- 
tion dimension Posy. Since Pos, is usually much smaller than 
Posy, we mainly focus on the effect of different Posy on the 
MEA pressure distribution, with Pos, being kept as constant. 
The translation directions of the four bolts position are shown 
by the arrows in Fig. 3. Because the end plate is foursquare and 
the translation of the four bolts are at the same time, the Posy of 
one bolt can represent the other three bolts’. 

In this study, the bolt position Pos, and the assembly pressure 
P were considered as randomly varying parameters with given 
distribution type and acted as the design variables. The normal 
stress and the normal stress uniformity were determined in terms 
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End Plate ; 
Translation Direction 


Fig. 3. Schematic of the translation of bolt position. 
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Table 2 

Probabilistic properties of the design variables and the coded design variables 
Variables* Lower Upper 
Pos, (mm) 0 33.5 

P (MPa) 0 1.5 
Coded Pos, (mm) —1.4285 1.4285 
Coded P (MPa) —1.4285 1.4285 


à All variables are assigned uniform distribution between the lower and upper 
limits, respectively. 


of the probabilistic design variables. The membrane’s maximum 
and minimum compression stress maxo; and ming, can be easily 
found. The difference between maximum and minimum com- 
pression stress can be defined as the differential compression 
stress in the membrane Ao. maxo; represents the value of con- 
tact pressure and Ao, represents the uniformity of MEA pressure 
distribution. When maxo; increases, it means the contact pres- 
sure on MEA becomes larger. If Ao, increases, the uniformity 
of MEA pressure distribution will become worse. 

According to the dimension of the end plate and the engineer- 
ing need, two design spaces for the two design variables were 
established, respectively. The random distribution types for the 
design variables must be chosen carefully in the design space 
so as to reveal the actual distributions of the design variables. 
Within a certain interval in the design space, we consider each 
value of the design variable makes the same contribution to the 
establishment of the response surface. That means any level of 
the design variables can be the robust solution equally. Since 
the “uniform distribution” is a very fundamental distribution for 
cases where there is no evidence that any value of the random 
variable is more likely than any other in the design space, we 
assign the uniform distribution for the design variables in the 
design space. Probabilistic properties of all the design variables 
are listed in Table 2. The response variables are maxo; and Aoz. 


4.2. Response surface establishment 


In order to find the relationship between the design variables 
and the response variables, the RSM was utilized. RSM is a 
combination of statistical and mathematical techniques, which is 
useful for developing, improving, and optimizing process [11]. 
A second-order model is commonly used for the multidisci- 
plinary design in RSM. In general, the response model can be 
written as follows: 


y(x) = Bo + x"b + x! Bx (1) 


where £o, b, B are coefficients or coefficient matrices for design 
variables x. 

The central composite design (CCD) is a popular method 
that can be effectively applied to construct the second-order 
model for RSM [12]. To fulfill the CCD, probabilistic design 
system (PDS) was used. PDS is an analysis technique for assess- 
ing the effect of uncertain input parameters and assumptions 
on the model. It can account for the randomness in input vari- 
ables such as material properties, boundary conditions, loads and 
geometry [13]. In PDS, statistical distribution functions (such as 


the Gaussian, the uniform distribution, etc.) describe uncertain 
parameters. For a given set of the distribution parameters of the 
design variables and the assumed distributions, one can generate 
a large set of random numbers for each variable using PDS [14]. 
In PDS, the CCD is fully automated and execution of CCD on 
the numerical model developed in Section 3 will result in a set of 
the response attributes. The results of CCD are listed in Table 3. 

Two second-order polynomial models in the form of Eq. (1) 
were fitted to the data. The items/effects that are not significant 
(confidence level=95%) were stepped down from the models 
using forward-stepwise-regression without damaging the model 
hierarchy. The values of regression coefficients were calculated 
and the fitted equations (in terms of coded values) for maxo, 
and Ao, are 


maxo; = 0.341 + 0.056x; + 0.236x2 + 0.043x1x2 (2) 
Ao, = 0.218 + 0.076x; + 0.158x2 + 0.016x7 + 0.057x1x2(3) 


where xı and x2 are the coded values of Pos, and P. The proba- 

bilistic properties of x; and x2 are listed in Table 2. The coding 

equations are 

xı = 0.085 Pos, — 1.429; x2 = 1.905 P — 1.429 (4) 
Analysis of variance (ANOVA) was conducted to determine 

the significance of the fitted regression models. Significance was 

judged by determining the probability level that the F-statistic 


calculated from the data is less than 5%. The F-statistic equation 
is 


__ SSR/(m — 1) 


~ SSE/(n — m) ©) 


where SSR is the sum of squares due to regression and SSE is the 
sum of squares of the residuals. n is the number of experiments 
and m is the number of terms in the fitted model. SSR and SSE 
are computed as 


SSR=5°G;-5);  SSE= S001 - 91)” (6) 
i=] i=l] 


where y; represents the ith response value predicted by the fitted 
model and y represents the average response value. 


Table 3 

Results of CCD 

Loop Posy (mm) P (MPa) maxo; (MPa) Ao, (MPa) 
1 16.75 0.7500 0.3305 0.2176 
2 0.1675 0.7500 0.2811 0.1497 
3 33.33 0.7500 0.4240 0.3481 
4 16.75 0.0075 0.0033 0.0022 
5 16.75 1.4925 0.6577 0.4330 
6 5.024 0.2250 0.0846 0.0460 
7 28.48 0.2250 0.1213 0.0947 
8 5.024 1.275 0.4797 0.2605 
9 28.48 1.275 0.6875 0.5364 


764 


The total variation of the response data, SST (total sum of 
squares) is computed as 


n 
SST = SSR + SSE = X `Q; - 5) (7) 


i=1 


R? is an accompanying statistic to the F-statistic. It expresses 
the proportion of the variation of y; from the model and the 
experimental data about the mean y. Another useful measure is 
the adjusted R? statistic Ri: 


Rg SSR, gp _, _ SSE/(n—m) 


~ SST’ A SST/(n— 1) (8) 


The ANOVA for the refined models is summarized in Table 4. 
Because of F(3, 5, 0.05) =5.41, F(4, 4, 0.05) = 6.39, the F-value 
of 796>F(3, 5, 0.05) for maxo, and 831>F(4, 4, 0.05) for 
Ao, imply the two models are significant at 95% confidence 
level. On the other hand, R? and RI were calculated to check 
the model adequacy. A high proportion of variability (R? > 0.9) 
in the response models can be explained successfully by the 
models (Eqs. (2) and (3)). However, a large value of R? does not 
always imply that the regression model is a good one. Adding 
a variable to the model will always increase R?, regardless of 
whether the additional variable is statistically significant or not. 
Thus, there is a need to use a Rx to evaluate the model adequacy 
and should be over 90%. Table 4 shows that R? and Rx values for 
the models do not differ dramatically. This can be an indication 
that insignificant terms have not been included in the models. 

To visualize the combined effect of the two design variables 
on the responses, the response surfaces and contour plots were 
generated for each of the fitted models. Fig. 4 shows the effect 
of Pos, and P on maxo,, and Ao;. 

As clear from Fig. 4(a), the maximum response for maxo, 
occurs when Posy and P are at their highest level. Increasing P 
rises up maxo; rapidly, which is in accordance with Lee’s results 
[6]. This suggests that the assembly pressure P has a very signif- 
icant effect on the MEA pressure distribution. It can be seen that 
the response also verifies noticeably at different levels of Posy, 
especially when the assembly pressure P is larger. With P fixed 
and Pos, increased, maxo; also becomes larger, which means the 
assembly bolts at the corner of end plate produce larger maxo; 
than other positions with the same assembly pressure. 

In this study, since Ao, represents the uniformity of the MEA 
pressure distribution and smaller value of Ao, means better uni- 
formity, the Ao, is subjected to be minimized. As shown in 
Fig. 4(b), the increase of Ao, is slower at the lower level of P 
and Posy than the higher level, and smaller value can be gained at 
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Fig. 4. Response surfaces and contour plots for (a) maxo; and (b) Ao;. 


the lower level of P and Posy. It can also be seen from Fig. 4 that 
the two response variables maxo; and Ao, are competing with 
each other, which means we cannot get a larger maxo; while 
keeping Ao, smaller simultaneously as we expect. 


4.3. Error propagation equations establishment for 
response variables 


The design may not achieve the desired result due to the influ- 
ence of uncontrollable noise and variation of the input factors 
[15]. Therefore, designs are sought that are not only optimal 
but also robust (insensitive) to inevitable changes in the noise 
and input factors. To accomplish this objective, one should set 
the controllable factors to the levels that reduce variation in the 
response: (1) caused by the variation in the uncontrollable fac- 
tors [16] and (2) transmitted from variation in the controllable 


Table 4 
ANOVA for the refined models 
Source maxo; Ao; 

d.f. Sum of squares F R? Ry d.f. Sum of squares F R? R: 
Regression 3 0.4772 4 0.2601 
Residual 5 9.99 x 1074 796 0.9979 0.9967 4 5.22 x 1074 831 0.9980 0.9968 
Total 8 0.4782 8 0.2606 
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factors [17]. In this research, we mainly focused on reducing 
transmitted variation in the controllable factors. 

Desirability-based robust design is a tool to find the control- 
lable factor settings that optimize the objective yet minimize the 
response variation of the design [15]. It requires construction of 
a response surface using a mathematical model (Eq. (1)). The 
variation transmitted to the response can be determined by the 
well-known error propagation equation, which is modeled by 
taking the partial derivatives of the polynomial (Eq. (1)) with 
respect to the design variables: 


1/2 


(9) 


n ay 2 
2 2 
by (>) Oxi + Oresid 


POE = oy = | 
i=1 


where oy is the model-predicted standard deviation of the 
response y, or, the variance of design variable x; and Oresid is 
the residual variance. 

The two POE for max o, and Ao, were calculated based 
on Eqs. (2), (3) and (9) with the standard deviation in P set as 
op =0.04 MPa, in Posy set as opos, = 0.33 mm and the residual 


variance set as 0: 


POE, = (3.253 + 1.171x1-+0.038x2+0.106x? + 0.015x2) 
x107? (10) 
POE, = (1.498 + 1.081x1 + 0.068x2 + 0.195x7 
+0.026x3 + 0.028x1x2) 7 x 1072 (11) 


where x; and x2 are the coded values of Posy and P. POE; and 
POE) are the standard deviation of maxo, and Ag, respecti- 
vely. 

To reduce the variance in the response, POE should be min- 
imized, therefore, it can be treated as an additional response 
built into the design process. The simultaneous optimization 
of several responses (in this case, maxo;z, Acz, POE, and 
POE?) is the essence of the desirability-based robust design 
[17,18]. 


4.4. Desirability function 


Desirability function is based on the idea that the “quality” 
of a product or process that has multiple quality characteristics, 
with one of them outside of some “desired” limits, is com- 
pletely unacceptable. The method finds operating conditions x 
that provide the “most desirable” response values. Depending 
on whether a particular response y; is to be maximized or min- 
imized, different desirability functions d;(y;) can be used. Let 
Li, T; and U; be the lower, target and upper values, respectively, 
that are desired for response y; with Li, T; and U; [19]. 

If a response is to be maximized, its individual desirability 
function is with the exponent S determining how important it is to 
hit the target value. For S= 1, the desirability function increases 
linearly towards T; which denotes a large enough value for the 
response; for S<1, the function is convex, and for S>1, the 


function is concave [19]: 


0, yi) < Li 
i(x) — Li\” 
di) = CRE) usas a2 
i yx) > T, 


If a response is to be minimized, its individual desirability func- 
tion is with 7; denoting a small enough value for the response: 


I yi) <T, 
i(x) — Ui \* 
di) = (eae). Tsss (13) 
0, yi(x) > Ui 


where T; represents a small enough value for the response. 

A single overall desirability index D can be obtained by com- 
bining the desirability value of each response. By assigning a 
range of numbers, for example 1-5, to the importance of opti- 
mizing each response variable, the weighing coefficients can be 
further refined. The final desirability index then is computed as 
follows: 


D = (d! x d? x d¥3 x- x duayl/ X wi 


n VOD Wi 
= (ie) (14) 
i=l 


where wj; is a number indicating the relative importance of the 
ith response, which might typically be an integer in the range of 
1-5, with 5 indicating the most importance and 1 indicating the 
least importance [19]. 

In this study, the two responses maxo; and Ao; were scaled 
to the desirability functions dı and d2 ranged from [0,1] by 
Eqs. (2), (3), (12) and (13) with Tı =0.4, Lı =0.8, T2 =0.1 and 
U2 =0.5, where maxo, and Ao; are subjected to be maximized 
and minimized, respectively. POE; and POE) were also scaled to 
two desirability functions d3 and d4 according to Eqs. (10), (11) 
and (13) with T3 = 0.01, U3 =0.02, T4 =0.005 and U4 =0.015. 

From Eq. (14), the overall desirability function for the robust 
design was then expressed as follows: 


Dyobust = (di x dy x d3 x dy)'/4 (15) 


where D,opust 18 the overall desirability function for the robust 
design. The weighing coefficients for all responses were selected 
as 1. 


5. Results and discussions 


Fig. 5(a) shows the overall desirability function Dyobust. It can 
be seen that there are a lot of zero values for some combines of 
P and Posy. This is because the responses of these combines are 
out of the “desired” limits of certain single desirability function 
d; (i= 1-4) and they are completely unacceptable for the overall 
desirability function. 

From Fig. 5, the overall desirability function is a continu- 
ous, nonlinear, piecewise function and therefore a direct search 
method, downhill simplex search, was used to find the combine 
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Table 5 
The statistic properties of maxo, and Ao; of three solutions 
P (MPa) Pos, (mm) maxo, (MPa) Ao; (MPa) 
Mean s* (x 1077) Mean sè (x 1077) 
Arbitrary solution 1.2 10 0.4893 1.63 0.2869 1 
Optimal solution 1.5 5.7759 0.5677 1.52 0.3104 0.87 
Robust solution 1.5 0 0.5098 1.38 0.2517 0.67 


è S.D. 


of the design variables Posy and P that maximizes the over- 
all desirability. By maximizing Drobust, the robust solution was 
obtained as — 1.4285 for coded value of Posy and 1.4285 for 
coded value of P, that is 0mm for Posy and 1.5 MPa for P with 
Dyobust = 9.5479 as Fig. 5(a) shows. 

To verify the reliability of the robust solution, we found the 
optimal solution of the assembly parameters, which will not con- 
sider the variations in the design variables. From the construction 
process of Drobust, one can see that the desirability function of 
the optimal design Dop for the criteria that maximum maxo; and 
minimum Ao, can be obtained by omitting d3 and d4, that is: 


Dop = (di x d2)? (16) 


where d and dz are the desirability functions of maxo, and øz. 


Desirability 


(a) Posy(mm) 


X: 5.776 
“YS 1.5 
Z: 0.4461 


Desirability 


(b) Posy(mm) 0 


P(MPa) 


Fig. 5. Contours of the overall desirability functions of (a) robust design and (b) 
optimal design. 


By maximizing Dop, the optimum solution was obtained as 
—0.9359 for coded value of Posy and 1.4285 for coded value of P, 
that is 5.776 mm for Posy and 1.5 MPa for P with Dop =0.4461. 
The overall desirability function of optimal solution Dop is 
shown in Fig. 5(b). 

An arbitrary combine of the design variables, Posy = 10 mm 
and P= 1.2 MPa, was also chosen to conduct the validation. 

In order to verify the robustness of the robust solution, we 
compared the performances of the arbitrary solution, optimal 
solution and robust solution on maxo, and Ao,. Three sets of 
normal distributions for P and Posy were established with the 
different mean values but the same variance. The mean values 
are the arbitrary solution, optimal solution and robust solu- 
tion, respectively, and the same variances are op = 0.04 MPa and 
oPosy = 0.33 mm. Using PDS again with the three normal dis- 
tributions of the design variables, three sets of maxo, and Ao, 
can be obtained corresponding to the arbitrary solution, optimal 
solution and robust solution, respectively. 

The statistic properties of maxo; and Ao, corresponding to 
the arbitrary solution, optimal solution and robust solution are 
listed in Table 5, respectively. As Table 5 shows, the robust 
solution shows less variability on the maxo, and Ao, than 
both the optimal solution and arbitrary solution. It can be also 
observed that the robustness has been achieved at the expense 
of the decrease of mean value of maxo; and Ao,. In this 
study, the larger maxo, and smaller Ao, the better. In this fash- 
ion, the robust solution decreases the value of contact pressure 
but makes the uniformity better. This can be explained by the 
competing relationship between the two responses maxo; and 
Aoz. 

Table 5 also shows the robustness of the robust solution, 
compared with the optimal solution, for Ag, is better than that 
for maxo,;. That is 10.14% standard deviation less on maxo, 
while 29.85% on Aoz. This is in accordance with Fig. 4, which 
also shows the gradient change of maxo, is less than that of 
Ao;. 


6. Conclusions 


A methodology of robust design integrated with FEA, RSM, 
POE and desirability function has been developed to find the 
robust solution of the assembly parameters for a PEMFC stack. 
It has been demonstrated that the proposed robust design scheme 
is feasible and effective for the assembly parameters consid- 
ered. For this PEMFC stack, the robust solution shows that the 
end plate bolts should be at the middle of four margins of the 
end plate and the assembly pressure should be 1.5 MPa. It also 
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shows that the assembly pressure plays a more important rule 
on the MEA pressure distribution than the position of end plate 
bolts. 

The RSM, PDS, POE and desirability function can be use- 
ful in the robust design process. Especially, the PDS, which 
integrated the probability and FEA, can be used to find out the 
influence of the variations of the design factors. The methodol- 
ogy of the present study can be used to seek robust sets of the 
assembly parameters for the designers to improve the perfor- 
mance of a PEMFC stack. 

In this study, the weighing coefficients for all response vari- 
ables were selected as the same values 1. However, from the 
construction process of the overall desirability function, it can 
be seen that one can change the weighing coefficients to vary 
the importance of each response variable individually based 
on individual engineering needs. The model was simplified by 
assuming isotropic material behavior, however, to make the 
results more accurate, the anisotropic material behavior should 
be researched and considered. The model also needs to be exper- 
imentally validated. Our research works in the future will be 
focused on these areas. 
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